L10: Geometric operations with rasters

Mosaicking, Resampling, Reprojecting

Bogdan G. Popescu

bogdan.popescu@johncabot.edu

John Cabot University

Mosicking Rasters

In this lectures, we will look at making changes to geometric component of rasters:

mosaicking
resampling
reprojecting

We will use packages like stars and units

Mosaicking: Elevation

In the next we few slides, we will use a Digitial Elevation Models (DEM) of the world, by mosaicking, subsetting, and reprojecting

We start putting all tiles together

Downloading DEM data

DEM stands for a digital elevation model

There are many data sources

The most helpful website for downloading DEM is: http://www.webgis.com/terr_world.html

Working with Elevation

Elevation

#Step1: Load the DEM files for every continent
#EUROPE
filepath<-"./data/elevation/w020n90/W020N90.DEM"
x1<-read_stars(filepath)
filepath<-"./data/elevation/w020n40/w020n40.DEM"
x2<-read_stars(filepath)
filepath<-"./data/elevation/e020n40/e020n40.DEM"
x3<-read_stars(filepath)
filepath<-"./data/elevation/e020n90/e020n90.DEM"
x4<-read_stars(filepath)
filepath<-"./data/elevation/w060n40/w060n40.DEM"
x5<-read_stars(filepath)
filepath<-"./data/elevation/w060n90/w060n90.DEM"
x6<-read_stars(filepath)

#ASIA
filepath<-"./data/elevation/e140n40/e140n40.DEM"
x7<-read_stars(filepath)
filepath<-"./data/elevation/e140n90/e140n90.DEM"
x8<-read_stars(filepath)
filepath<-"./data/elevation/e060n40/e060n40.DEM"
x9<-read_stars(filepath)
filepath<-"./data/elevation/e060n90/e060n90.DEM"
x10<-read_stars(filepath)
filepath<-"./data/elevation/e100n40/e100n40.DEM"
x11<-read_stars(filepath)
filepath<-"./data/elevation/e100n90/e100n90.DEM"
x12<-read_stars(filepath)

#AUSTRALIA
filepath<-"./data/elevation/e060s10/e060s10.DEM"
x13<-read_stars(filepath)
filepath<-"./data/elevation/e100s10/e100s10.DEM"
x14<-read_stars(filepath)
filepath<-"./data/elevation/e140s10/e140s10.DEM"
x15<-read_stars(filepath)
filepath<-"./data/elevation/e060s60/e060s60.DEM"
x16<-read_stars(filepath)
filepath<-"./data/elevation/e120s60/e120s60.DEM"
x17<-read_stars(filepath)

#NORTH AMERICA
filepath<-"./data/elevation/w180n90/w180n90.DEM"
x18<-read_stars(filepath)
filepath<-"./data/elevation/w140n90/w140n90.DEM"
x19<-read_stars(filepath)
filepath<-"./data/elevation/w100n90/w100n90.DEM"
x20<-read_stars(filepath)
filepath<-"./data/elevation/w180n40/w180n40.DEM"
x21<-read_stars(filepath)
filepath<-"./data/elevation/w140n40/w140n40.DEM"
x22<-read_stars(filepath)
filepath<-"./data/elevation/w100n40/w100n40.DEM"
x23<-read_stars(filepath)

#SOUTH AMERICA
filepath<-"./data/elevation/w180s10/w180s10.DEM"
x24<-read_stars(filepath)
filepath<-"./data/elevation/w140s10/w140s10.DEM"
x25<-read_stars(filepath)
filepath<-"./data/elevation/w100s10/w100s10.DEM"
x26<-read_stars(filepath)
filepath<-"./data/elevation/w180s60/w180s60.DEM"
x27<-read_stars(filepath)
filepath<-"./data/elevation/w120s60/w120s60.DEM"
x28<-read_stars(filepath)

#SOUTH AFRICA
filepath<-"./data/elevation/w060s10/w060s10.DEM"
x29<-read_stars(filepath)
filepath<-"./data/elevation/w020s10/w020s10.DEM"
x30<-read_stars(filepath)
filepath<-"./data/elevation/e020s10/e020s10.DEM"
x31<-read_stars(filepath)
filepath<-"./data/elevation/w060s60/w060s60.DEM"
x32<-read_stars(filepath)
filepath<-"./data/elevation/w000s60/w000s60.DEM"
x33<-read_stars(filepath)

Elevation

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Lower dx values. Lower dx values means lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
st_crs(elev)<-st_crs(borders)

fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Digital Elevation Model (DEM)

Digital Elevation Model (DEM) in Europe

library(purrr)

#Identifying coordinates of Northern points
norway<-subset(borders, SOVEREIGNT == "Norway")
nc_geom <- st_geometry(norway)
list_points_norway<-nc_geom[[1]]
list_points_norway<-flatten(list_points_norway)
list_points_norway2<-as.data.frame(do.call(rbind, list_points_norway))
names(list_points_norway2)<-c("lon_x", "lat_y")
max_lat_y_eu<-max(list_points_norway2$lat_y)

#Identifying coordinates of Southern points
greece<-subset(borders, SOVEREIGNT == "Greece")
nc_geom <- st_geometry(greece)
list_points_greece<-nc_geom[[1]]
list_points_greece<-flatten(list_points_greece)
list_points_greece2<-as.data.frame(do.call(rbind, list_points_greece))
names(list_points_greece2)<-c("lon_x", "lat_y")
min_lat_y_eu<-min(list_points_greece2$lat_y)

#Identifying coordinates of Eastern points
ukraine<-subset(borders, SOVEREIGNT == "Ukraine")
nc_geom <- st_geometry(ukraine)
list_points_ukraine<-nc_geom[[1]]
list_points_ukraine<-flatten(list_points_ukraine)
list_points_ukraine2<-as.data.frame(do.call(rbind, list_points_ukraine))
names(list_points_ukraine2)<-c("lon_x", "lat_y")
max_lon_x_eu<-max(list_points_ukraine2$lon_x)

#Identifying coordinates of Western points
portugal<-subset(borders, SOVEREIGNT == "Portugal")
nc_geom <- st_geometry(portugal)
list_points_portugal<-nc_geom[[1]]
list_points_portugal<-flatten(list_points_portugal)
list_points_portugal2<-as.data.frame(do.call(rbind, list_points_portugal))
names(list_points_portugal2)<-c("lon_x", "lat_y")
min_lon_x_eu<-min(list_points_portugal2$lon_x)

Digital Elevation Model (DEM) in Europe

Digital Elevation Model (DEM) in the US

Elevation 3D

Raster Resampling

Raster resampling is the process of transferring raster values from the original grid to a different grid.

Resampling may be necessary for:

alligning several input rasters that come from different sources to the same grid
Reducing the resolution of very detailed rasters, so that they are more convenient to work with

Raster Resampling

In case you have not noticed, we have already resampled the DEM raster.

#Step2: Creating a mosaic
star_mosaic_eu <- st_mosaic(x1, x2, x3, x4, x5, x6)
star_mosaic_rus <- st_mosaic(x7, x8, x9, x10, x11, x12)
star_mosaic_aus <- st_mosaic(x13, x14, x15, x16, x17)
star_mosaic_namer <- st_mosaic(x18, x19, x20, x21, x22, x23)
star_mosaic_samer <- st_mosaic(x24, x25, x26, x27, x28)
star_mosaic_safr <- st_mosaic(x29, x30, x31, x32, x33)

#Step3: Creating a second mosaic
all<-st_mosaic(star_mosaic_eu, star_mosaic_rus, star_mosaic_aus, star_mosaic_namer, star_mosaic_samer, star_mosaic_safr)

#Step4: Lower dx values. Lower dx values means lower resolution
elev<-st_warp(all, st_as_stars(st_bbox(all), dx = 0.4))

#Step5: Renaming the raster file
names(elev)<-"elev"

#Step6: Ensuring the the raster has the same crs
st_crs(elev)<-st_crs(borders)

fig<-ggplot()+
  geom_stars(data=elev)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())

Nearest Neighbor Resampling

There are a variety of other resampling methods

In the nearest neighbor sampling, the new raster pixels get the value from the nearest pixel of the original raster.

Some values may be lost.

# Nearest Neighbor Resampling

#Step4: Lower dx values. Lower dx values means lower resolution
grid = st_as_stars(st_bbox(all), dx = 0.4, dy = 0.4)
dem1 = st_warp(elev, grid, method = "near")

fig<-ggplot()+
  geom_stars(data=dem1)+
  scale_fill_gradientn(colours=c("chartreuse3","red"),   na.value = NA)+
  geom_sf(data=borders, linewidth = 0.01, fill = NA, color = "white", alpha=0.5)+
  theme_bw()+
  theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank())
fig

Bilinear resampling

In bilinear resampling, each new raster cell gets a weighted average of four nearest cells from the input, rather than just one.

Bilinear resampling

Two technical things to note about st_warp when using a method other than method="near":

use_gdal=TRUE needs to be specified, otherwise the method argument is ignored
no_data_value needs to be set to a value outside of the valid range (such as no_data_value=-9999 for a DEM raster);

dem2 = st_warp(elev, grid, method = "bilinear", use_gdal = TRUE, no_data_value = -9999)

Average resampling

Another useful method is the average resampling method, where each new cell gets the weighted average of all overlapping input cells:

Average resampling

dem3 = st_warp(elev, grid, method = "average", use_gdal = TRUE, no_data_value = -9999)